
path1='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\Registration and analysis\4. position measurements\';
dir_to_search=[path1];
pathtemp='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\Registration and analysis\3.Channel colocalization\s1b2\';
path2='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\Registration and analysis\2.registered after downscaling to half res\';
path3='C:\Users\Tom�s\Documents\1. Regehr Lab\Data\HCR\'

% get the folder contents
d = dir(path1);
% remove all files (isdir property is 0)
dfolders = d([d(:).isdir]==1);
% remove '.' and '..'
dfolders = dfolders(~ismember({dfolders(:).name},{'.','..'}));

clear Cells ch1values ch2values ch3values

for j=1:numel(dfolders)

    
    
    %use vertcat 
    
    %small test for one image to see how well this works

    Cells=HCR(j).MLI2(:,[1:2]);

% Cells=HCR(1).Lugaro(:,[1:2]);
% ch1=imread([pathtemp 's1b2_CCs composite.tif'],1);
% ch2=imread([pathtemp 's1b2_CCs composite.tif'],2);
% ch3=imread([pathtemp 's1b2_CCs composite.tif'],3);

ch1=imread([path2 dfolders(j).name '\' 'registered_' dfolders(j).name '_Nxph1.tiff'],1);
ch2=imread([path2 dfolders(j).name '\' 'registered_' dfolders(j).name '_Slc6a5.tiff'],1);
ch3=imread([path2 dfolders(j).name '\' 'registered_' dfolders(j).name '_Aldh1a3.tiff'],1);



% figure
% imagesc(ch1)
% caxis([0 300])
% hold on
% scatter(CCs(:,1),CCs(:,2),10,'filled','r')

ch1Mat=nan(61,61,size(Cells,1));
ch2Mat=nan(61,61,size(Cells,1));
ch3Mat=nan(61,61,size(Cells,1));
for i=2:size(Cells,1)
    ch1Mat(:,:,i)=ch1([round(Cells(i,2))-30:round(Cells(i,2))+30],round([round(Cells(i,1))-30:round(Cells(i,1))+30]));
    ch2Mat(:,:,i)=ch2([round(Cells(i,2))-30:round(Cells(i,2))+30],round([round(Cells(i,1))-30:round(Cells(i,1))+30]));
    ch3Mat(:,:,i)=ch3([round(Cells(i,2))-30:round(Cells(i,2))+30],round([round(Cells(i,1))-30:round(Cells(i,1))+30]));

end

if j==1
    cumch1Mat=ch1Mat;
    cumch2Mat=ch2Mat;
    cumch3Mat=ch3Mat;
else
end

% concatenate across slices
if ~(j==1)
    cumch1Mat=cat(3,cumch1Mat,ch1Mat);
    cumch2Mat=cat(3,cumch2Mat,ch2Mat);
    cumch3Mat=cat(3,cumch3Mat,ch3Mat);
else
end





end


% %% Display  mosaic of top cells (3 rows, 3 colums, total 9 cells)
% for idx=1:9
% randidx(idx)=round(rand(1)*size(cumch1Mat,3));
% end
% 
% 
% 
% row1=cat(1,cat(1,cumch1Mat(:,:,randidx(1)),cumch1Mat(:,:,randidx(2)),cumch1Mat(:,:,randidx(3))));
% row2=cat(1,cat(1,cumch1Mat(:,:,randidx(4)),cumch1Mat(:,:,randidx(5)),cumch1Mat(:,:,randidx(6))));
% % row3=cat(1,cat(1,cumch1Mat(:,:,randidx(7)),cumch1Mat(:,:,randidx(8)),cumch1Mat(:,:,randidx(9))));
% 
% mosaic1=cat(2,row1,row2);
% 
% row1=cat(1,cat(1,cumch2Mat(:,:,randidx(1)),cumch2Mat(:,:,randidx(2)),cumch2Mat(:,:,randidx(3))));
% row2=cat(1,cat(1,cumch2Mat(:,:,randidx(4)),cumch2Mat(:,:,randidx(5)),cumch2Mat(:,:,randidx(6))));
% % row3=cat(1,cat(1,cumch2Mat(:,:,randidx(7)),cumch2Mat(:,:,randidx(8)),cumch2Mat(:,:,randidx(9))));
% 
% mosaic2=cat(2,row1,row2);
% 
% row1=cat(1,cat(1,cumch3Mat(:,:,randidx(1)),cumch3Mat(:,:,randidx(2)),cumch3Mat(:,:,randidx(3))));
% row2=cat(1,cat(1,cumch3Mat(:,:,randidx(4)),cumch3Mat(:,:,randidx(5)),cumch3Mat(:,:,randidx(6))));
% % row3=cat(1,cat(1,cumch3Mat(:,:,randidx(7)),cumch3Mat(:,:,randidx(8)),cumch3Mat(:,:,randidx(9))));
% 
% mosaic3=cat(2,row1,row2);
% 
% % 
% figure;
% imagesc(mosaic1)
% colormap([cmap('green',100,0,50)]);
% caxis([200 750])
% 
% figure
% imagesc(mosaic2)
% colormap([cmap('cyan',100,0,50)]);
% caxis([200 1500])
% 
% figure
% imagesc(mosaic3)
% colormap([cmap('magenta',100,0,50)]);
% caxis([200 1500])

% 
 %%
% imwrite(uint16(mosaic1),[path3 'mosaic.tif'])
% imwrite(uint16(mosaic2),[path3 'mosaic.tif'],'WriteMode','append')
% imwrite(uint16(mosaic3),[path3 'mosaic.tif'],'WriteMode','append')

%%
% make nice figure 
%display all the cells stacked together to look at data distribution
radius=16;
ylimitprofiles=370;
figure
c1Ax=subaxis(2,3,1)
  imagesc(nanmean(cumch1Mat,3))
  caxis([0 ylimitprofiles/2])
    y1=colormap(c1Ax,[cmap('green',100,0,50)]);

  axis square
  axis off
  subaxis(2,3,4)
%   plot(mean(mean(cumch1Mat,3),1))
meanGreen=nanmean(cumch1Mat,3);
greenHor=mean(meanGreen([29:33],:),1);
greenVert=mean(meanGreen(:,[29:33]),2);
greenProf=mean([greenVert';greenHor],1)
  hold on
  plot(greenProf,'k')
%   plot(mean(mean(cumch1Mat,3),2)) %average profiles in x and y to decide on boundary of circular mask
  
  ylim([0 ylimitprofiles/2])
%   vline([30-radius 30+radius],'--k')
  %draw axis for visual purpose
  vline(0,'k')
  hline(0,'k')
  axis square
  axis off
  
 c2Ax= subaxis(2,3,2)
  imagesc(nanmean(cumch2Mat,3))
  caxis([0 ylimitprofiles])
  y2=colormap(c2Ax,[cmap('cyan',100,0,50)]);
  axis square
  axis off
  subaxis(2,3,5)
%   plot(mean(mean(cumch2Mat,3),1))
%   hold on
%   plot(mean(mean(cumch2Mat,3),2)) %average profiles in x and y to decide on boundary of circular mask
meanCyan=nanmean(cumch2Mat,3);
CyanHor=mean(meanCyan([29:33],:),1);
CyanVert=mean(meanCyan(:,[29:33]),2);
CyanProf=mean([CyanVert';CyanHor],1)
  hold on
  plot(CyanProf,'k')
   ylim([0 ylimitprofiles])
%    vline([30-radius 30+radius],'--k')
   %draw axis for visual purpose
  vline(0,'k')
  hline(0,'k')
  axis square
  axis off
  
   c3Ax= subaxis(2,3,3)
  imagesc(mean(cumch3Mat,3))
  
  caxis([0 ylimitprofiles])
    y3=colormap(c3Ax,[cmap('magenta',100,0,50)]);

  axis square
  axis off
  subaxis(2,3,6)
%   plot(mean(mean(cumch3Mat,3),1))
%   hold on
%   plot(mean(mean(cumch3Mat,3),2)) %average profiles in x and y to decide on boundary of circular mask
meanMagenta=nanmean(cumch3Mat,3);
MagentaHor=mean(meanMagenta([29:33],:),1);
MagentaVert=mean(meanMagenta(:,[29:33]),2);
MagentaProf=mean([MagentaVert';MagentaHor],1)
  hold on
  plot(MagentaProf,'k')
   ylim([0 ylimitprofiles])
%    vline([30-radius 30+radius],'--k')
   %draw axis for visual purpose
  vline(0,'k')
  hline(0,'k')
  axis square
  axis off
 
  %% masking
    
% Create a logical image of a circle with specified
% diameter, center, and image size.
% First create the image.
imageSizeX = 61;
imageSizeY = 61;
[columnsInImage, rowsInImage] = meshgrid(1:imageSizeX, 1:imageSizeY);
% Next create the circle in the image.
centerX = ceil(imageSizeX / 2); % in the middle.
centerY = ceil(imageSizeY / 2);
radius = 25; %based on profiles

circlePixelsLogical = (rowsInImage - centerY).^2 ...
    + (columnsInImage - centerX).^2 <= radius.^2;
% circlePixels is a 2D "logical" array.
% Now, display it.
figure
image(circlePixelsLogical) ;
colormap([0 0 0; 1 1 1]);
title('Binary image of a circle');
axis square;

circlePixels01=double(circlePixelsLogical);
test=circlePixels01.*cumch1Mat; %mask out crap outisde cells 
figure
imagesc(nanmean(test,3))
caxis([0 50])
axis square

%% take averages of each cell for each channel inside the circular mask
clear ch1values
clear ch2values
clear ch3values
clear ch1tsne
clear ch2tsne
clear ch3tsne
for c=1:size(cumch1Mat,3)
    tempdata1=cumch1Mat(:,:,c);
%     ch1tsne(c,:)=tempdata1(circlePixelsLogical);
    ch1values(c)=mean(tempdata1(circlePixelsLogical));
    
     tempdata2=cumch2Mat(:,:,c);
%       ch2tsne(c,:)=tempdata2(circlePixelsLogical);
    ch2values(c)=mean(tempdata2(circlePixelsLogical));
    
     tempdata3=cumch3Mat(:,:,c);
%       ch3tsne(c,:)=tempdata3(circlePixelsLogical);
    ch3values(c)=mean(tempdata3(circlePixelsLogical));
    
end

figure
limVal=300;
subaxis(1,3,1)
scatter(ch1values,ch2values,10,'filled','k')
xlabel('Nxph1')
ylabel('Slc6a5')
axis square
xlim([0 limVal])
ylim([0 limVal])
subaxis(1,3,2)
scatter(ch1values,ch3values,10,'filled','k')
xlabel('Nxph1')
ylabel('Aldh1a3')
xlim([0 limVal])
ylim([0 limVal])
axis square
subaxis(1,3,3)
scatter(ch2values,ch3values,10,'filled','k')
xlabel('Slc6a5')
ylabel('Aldh1a3')
xlim([0 limVal])
ylim([0 limVal])
axis square




%% manually add values to cell types for cumulative plot (this is lazy but only needs tob e done once).

% 
% CCValues(:,1)=ch1values;
% CCValues(:,2)=ch2values;
% CCValues(:,3)=ch3values;
% 
% 
% GlobularValues(:,1)=ch1values;
% GlobularValues(:,2)=ch2values;
% GlobularValues(:,3)=ch3values;

% 
% 
% LugaroValues(:,1)=ch1values;
% LugaroValues(:,2)=ch2values;
% LugaroValues(:,3)=ch3values;


% 
% GolgiValues(:,1)=ch1values;
% GolgiValues(:,2)=ch2values;
% GolgiValues(:,3)=ch3values;


% MLI2Values(:,1)=ch1values;
% MLI2Values(:,2)=ch2values;
% MLI2Values(:,3)=ch3values;
% 


%% Make cumulative plots
MLI2Color=[233 83 30]./255; % MLI2 color
CCColor=[30, 180, 233]./255; % CCs color from illustrator
GlobularColor=[74 184 89]./255; %Globular color

LugaroColor=[127 63 152]./255; %from Illustrator
 GolgiColor=[177 179 182]./255; %Golgi color illustrator



binWidth=2;
figure
hold on
Ch1CChist=histogram(CCValues(:,1),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
Ch1Globularhist=histogram(GlobularValues(:,1),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
Ch1Lugarohist=histogram(LugaroValues(:,1),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
Ch1Golgihist=histogram(GolgiValues(:,1),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
Ch1MLI2hist=histogram(MLI2Values(:,1),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% Globularhist=histogram(GlobularPCLdistance./GlobularMLlength,'binwidth',binWidth,'FaceAlpha',0.5,'LineStyle', 'none')



figure
hold on
Ch1CCStair=stairs(Ch1CChist.BinEdges(1:end-1), Ch1CChist.Values,'color',CCColor,'linewidth',2)
Ch1GlobularStair=stairs(Ch1Globularhist.BinEdges(1:end-1), Ch1Globularhist.Values,'color',GlobularColor,'linewidth',2)
Ch1LugaroStair=stairs(Ch1Lugarohist.BinEdges(1:end-1), Ch1Lugarohist.Values,'color',LugaroColor,'linewidth',2)
Ch1GolgiStair=stairs(Ch1Golgihist.BinEdges(1:end-1), Ch1Golgihist.Values,'color',GolgiColor,'linewidth',2)
Ch1MLI2Stair=stairs(Ch1MLI2hist.BinEdges(1:end-1), Ch1MLI2hist.Values,'color',MLI2Color,'linewidth',2)


figure
hold on
Ch2CChist=histogram(CCValues(:,2),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
Ch2Globularhist=histogram(GlobularValues(:,2),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
Ch2Lugarohist=histogram(LugaroValues(:,2),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
Ch2Golgihist=histogram(GolgiValues(:,2),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
Ch2MLI2hist=histogram(MLI2Values(:,2),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% Globularhist=histogram(GlobularPCLdistance./GlobularMLlength,'binwidth',binWidth,'FaceAlpha',0.5,'LineStyle', 'none')



figure
hold on
Ch2CCStair=stairs(Ch2CChist.BinEdges(1:end-1), Ch2CChist.Values,'color',CCColor,'linewidth',2)
Ch2GlobularStair=stairs(Ch2Globularhist.BinEdges(1:end-1), Ch2Globularhist.Values,'color',GlobularColor,'linewidth',2)
Ch2LugaroStair=stairs(Ch2Lugarohist.BinEdges(1:end-1), Ch2Lugarohist.Values,'color',LugaroColor,'linewidth',2)
Ch2GolgiStair=stairs(Ch2Golgihist.BinEdges(1:end-1), Ch2Golgihist.Values,'color',GolgiColor,'linewidth',2)
Ch2MLI2Stair=stairs(Ch2MLI2hist.BinEdges(1:end-1), Ch2MLI2hist.Values,'color',MLI2Color,'linewidth',2)


figure
hold on
Ch3CChist=histogram(CCValues(:,3),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
Ch3Globularhist=histogram(GlobularValues(:,3),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
Ch3Lugarohist=histogram(LugaroValues(:,3),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
Ch3Golgihist=histogram(GolgiValues(:,3),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
Ch3MLI2hist=histogram(MLI2Values(:,3),'binwidth',binWidth,'FaceAlpha',0,'LineStyle', 'none','normalization' ,'cdf')
% Globularhist=histogram(GlobularPCLdistance./GlobularMLlength,'binwidth',binWidth,'FaceAlpha',0.5,'LineStyle', 'none')



figure
hold on
Ch3CCStair=stairs(Ch3CChist.BinEdges(1:end-1), Ch3CChist.Values,'color',CCColor,'linewidth',2)
Ch3GlobularStair=stairs(Ch3Globularhist.BinEdges(1:end-1), Ch3Globularhist.Values,'color',GlobularColor,'linewidth',2)
Ch3LugaroStair=stairs(Ch3Lugarohist.BinEdges(1:end-1), Ch3Lugarohist.Values,'color',LugaroColor,'linewidth',2)
Ch3GolgiStair=stairs(Ch3Golgihist.BinEdges(1:end-1), Ch3Golgihist.Values,'color',GolgiColor,'linewidth',2)
Ch3MLI2Stair=stairs(Ch3MLI2hist.BinEdges(1:end-1), Ch3MLI2hist.Values,'color',MLI2Color,'linewidth',2)



















% CCch1=ch1tsne; 156
% CCch2=ch2tsne; 55
% CCch3=ch3tsne; 607

% MLI2tsne=horzcat(ch1tsne./156,ch2tsne./55,ch3tsne./607);
% Golgitsne=horzcat(ch1tsne./156,ch2tsne./55,ch3tsne./607);
% CCtnse=horzcat(CCch1./156,CCch2./55,CCch3./607);
% Globulartsne=horzcat(ch1tsne./156,ch2tsne./55,ch3tsne./607); %uses masked circle pixels
%% Use t-sne on the masked images with appropriate labels too see clustering of different cell types.
% alternatively use average value of masked image for each channel and
% normalized position in ML.

%make X vector for tsne, it will be size(X)= numberofCells by linearized
%pixelspercell
%group matrix g is going to be size = number of cells by 1
% X=vertcat(CCtsne,Globulartsne,Golgitsne,MLI2tsne);
% g=vertcat(ones(size(CCtsne,1),1),2*ones(size(Globulartsne,1),1),3*ones(size(Golgitsne,1),1),4*ones(size(MLI2tsne,1),1))
% rng default % for reproducibility
% Y = tsne(X,'Algorithm','barneshut','NumPCAComponents',0,'perplexity',500);
% figure;
% gscatter(Y(:,1),Y(:,2),g)
%  
    